#Nature's Kidneys: the Role of Wetland Reserve Easements in Restoring Water Quality
#Nicole Karwowski and Marin Skidmore
#Replication File
#2025

#Clean and Merge the NRCS easement data to Subwatersheds

#open libraries
library(fst)
library(ggplot2)
library(raster)
library(sf)
library(dplyr)
library(readr)
library(mapview)
library(haven)
library(stringr)

########################################################################################

#set wd
setwd("C:/Users/16308/OneDrive - Montana State University/UW_Madison/Water quality/Replication Package")

#Clean WRP detailed easement data
#open data
Ease <- read_dta("Data/NRCS easements/Raw/NRCSEasementsParcels_20210914.dta")

#clean out the character NAs with NA
Ease$RestorationFY<-ifelse(Ease$RestorationFY=="NA", NA, Ease$RestorationFY)
Ease$RecordedFY<-ifelse(Ease$RecordedFY=="NA", NA, Ease$RecordedFY)

Ease$NationalInitiativeGroup<-ifelse(Ease$NationalInitiativeGroup=="NA", NA, Ease$NationalInitiativeGroup)
table(Ease$NationalInitiativeGroup, useNA="always")

#Dates
Ease$RestorationCompletionDate<-as.Date(Ease$RestorationCompletionDate, origin="1960-01-01")
Ease$EasementRecordedDate<-as.Date(Ease$EasementRecordedDate, origin="1960-01-01")

Ease$ClosingDate<-as.Date(Ease$ClosingDate, origin="1960-01-01")
Ease$ClosingFY<-substring(Ease$ClosingDate, 1,4)

Ease$ApplicationDate<-as.Date(Ease$ApplicationDate, origin="1960-01-01")
Ease$ApplicationFY<-substring(Ease$ApplicationDate, 1, 4)

Ease$EasementRecordedDate<-as.Date(Ease$EasementRecordedDate, origin="1960-01-01")
Ease$RecordedFY<-substring(Ease$EasementRecordedDate,1,4)

Ease$AgreementNumber<-str_trim(Ease$AgreementNumber)

Ease$AgreementStartFY<-Ease$AgeementStartFY


##YearMonthVariables
Ease$ApplicationYM<-substring(Ease$ApplicationDate, 1, 7)
Ease$ApplicationYM<-str_replace_all(Ease$ApplicationYM, "-", "")

Ease$ClosingYM<-substring(Ease$ClosingDate, 1, 7)
Ease$ClosingYM<-str_replace_all(Ease$ClosingYM, "-", "")

Ease$RecordedYM<-substring(Ease$EasementRecordedDate, 1, 7)
Ease$RecordedYM<-str_replace_all(Ease$RecordedYM, "-", "")

Ease$RestorationDoneYM<-substring(Ease$RestorationCompletionDate, 1, 7)
Ease$RestorationDoneYM<-str_replace_all(Ease$RestorationDoneYM, "-", "")

Details<-subset(Ease, select=c(Program, FPE, WR, State, CountyName, HUC, StateAbv, Fips, 
                               AgreementNumber, EnrollmentType,
                               ThirtyYear, TenYear, Permanent, CostShare, 
                               Duration, AgreementStatus, AgreementParcelAcres, ClosedAcres, 
                               AgreementParcelOrClosedAcres, AgreePercentageCrop,
                               AgreementParcelCropAcres, FinalAcres,
                               ApplicationFY, AgreementStartFY, EnrollmentFY,
                               ClosingFY, RecordedFY, RestorationFY,
                               ApplicationDate, ClosingDate, EasementRecordedDate,  
                               RestorationCompletionDate, ApplicationYM, ClosingYM,
                               RecordedYM, RestorationDoneYM, 
                               Latitude, Longitude))
#sort by date
Details<-Details[order(Details$ClosingDate, Details$FinalAcres),]

#create a row id
Details<-mutate(Details, Id=row_number())

Details<-Details %>% relocate(Id, .before=Program)

#fix the cropland percentage variable
Details$AgreePercentageCrop<-ifelse(Details$AgreePercentageCrop>1, 
                                       1, Details$AgreePercentageCrop)

Details$AgreementParcelCropAcres<-Details$AgreePercentageCrop*Details$AgreementParcelAcres

Details$FinalCropAcres<-Details$AgreePercentageCrop*Details$FinalAcres
Details<- Details %>% relocate(FinalCropAcres, .after=FinalAcres)

#take out the ten year contracts
Details<-subset(Details, Details$TenYear==0)

#take out the thirty year contracts
#use only permanenet easements in the analyses
Details<-subset(Details, Details$EnrollmentType=="Easement")

#impute the 2 missing application dates
Details$diff<-as.numeric(Details$ClosingDate-Details$ApplicationDate)
#summary(Details$diff)
#on average it takes 778 days from application to closing (26 months)
Details$ApplicationYM<-ifelse(Details$Id=="1919", "199512", Details$ApplicationYM)
Details$ApplicationYM<-ifelse(Details$Id=="3794", "199804", Details$ApplicationYM)

#impute the missing restoration dates
#on average, it takes 913 days between closing and restoration. 
#that is equal to 30 months, or 2.5 years. 
MissingDate<-subset(Details, is.na(RestorationCompletionDate))

MissingDate$RestorationCompletionDate<-MissingDate$ClosingDate + 913

Details<-subset(Details, !is.na(RestorationCompletionDate))

Details<-rbind(Details, MissingDate)

Details<-subset(Details, select=-c(diff))

#lets save this data set
save(Details, file="Data/NRCS easements/Processed/EasementsDetailedCleaned.RData")

############################################################################################

#open the huc12 data to create huc12 month year panel of WRP
huc12<-st_read("Data/WatershedBoundaries/Raw/USGS_WBD/USGS_WBD.shp")
huc12<-huc12[, c(6,8.20)]

#turn Details into SF with lat and long info 
DetailsHuc<-st_as_sf(Details, coords=c("Longitude", "Latitude"), crs=crs(huc12))

#check lat longs if they look right
mapview(DetailsHuc)

#join easements with huc 12 identifier
DetailsHuc<-st_join(DetailsHuc, huc12, join=st_within)

DetailsHuc<-DetailsHuc %>% rename(huc8=HUC)

DetailsHuc<-subset(as.data.frame(DetailsHuc), select=-geometry)

#Aggregate by huc12 year and month

#closing year month agg
ClosedEase<-aggregate(FinalAcres~huc12+ClosingYM, data=DetailsHuc, na.rm=TRUE, FUN=sum)
ClosedCropEase<-aggregate(FinalCropAcres~huc12+ClosingYM, data=DetailsHuc, na.rm=TRUE, FUN=sum)

closed<-cbind(ClosedEase, ClosedCropEase)
closed<-closed[,c(1:3, 6)]
names(closed)<-c("huc12", "YearMonth", "ClosedEasementAcres", "ClosedEasementCropAcres")

table(closed$ClosedEasementAcres>=closed$ClosedEasementCropAcres)

#restoration done year month agg
RestoredEase<-aggregate(FinalAcres~huc12+RestorationDoneYM, data=DetailsHuc, na.rm=TRUE, FUN=sum)
RestoredCropEase<-aggregate(FinalCropAcres~huc12+RestorationDoneYM, data=DetailsHuc, na.rm=TRUE, FUN=sum)

DetailsHuc$Project<-1
RestoredProjYM<-aggregate(Project~huc12+RestorationDoneYM, data=DetailsHuc, na.rm=TRUE, FUN=sum)

restored<-cbind(RestoredEase, RestoredCropEase)
restored<-cbind(restored, RestoredProjYM)
restored<-restored[,c(1:3,6,9)]

table(restored$FinalAcres>=restored$FinalCropAcres)

names(restored)<-c("huc12", "YearMonth", "RestoredEasementAcres", "RestoredEasementCropAcres", "RestoredProjects")

#application done 
ApplyEase<-aggregate(FinalAcres~huc12+ApplicationYM, data=DetailsHuc, na.rm=TRUE, FUN=sum)
ApplyCropEase<-aggregate(FinalCropAcres~huc12+ApplicationYM, data=DetailsHuc, na.rm=TRUE, FUN=sum)

applied<-cbind(ApplyEase, ApplyCropEase)
applied<-applied[,c(1:3,6)]

names(applied)<-c("huc12", "YearMonth", "ApplyEasementAcres", "ApplyEasementCropAcres")

#merge them together
EaseAll<-merge(applied, closed, by=c("huc12", "YearMonth"), all.x=TRUE, all.y=TRUE)
EaseAll<-merge(EaseAll, restored, by=c("huc12", "YearMonth"), all.x=TRUE, all.y=TRUE)

#make the NAs zero
EaseAll[is.na(EaseAll)]<-0

#this represents the new easements in each huc 12 year month
ease_huc12<-EaseAll

save(ease_huc12, file="Data/NRCS easements/Processed/easements_panel_huc_ym.RData")



